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We study the dynamics of a collection of nonlinearly coupled limit cycle oscillators relevant to a wide class of 
systems, ranging from neuronal populations to electrical circuits, over network topologies varying from a 
regular ring to a random network. We find that for sufficiently strong coupling strengths the trajectories of 
the system escape to infinity in the regular ring network. However when a fraction of the regular connections 
are dynamically randomized, the unbounded growth is suppressed and the system remains bounded. 
Further, we find a scaling relation between the critical fraction of random links necessary for successful 
prevention of explosive behavior and the network rewiring time-scale. These results suggest a mechanism by 
which blow-ups may be controlled in extended oscillator systems. 

Collections of coupled dynamical systems were introduced as simple models capturing the essential features 
of large interactive systems'. Research focused on such distributed systems and complex networks has 
provided frameworks for understanding spatiotemporal patterns in problems ranging from Josephson 
junction arrays and multimode lasers, to microfluidic arrays and evolutionary biology'"^. 

In this work we study the dynamics of nonlinearly coupled oscillators relevant to systems ranging from 
neuronal populations to electrical circuits, over network topologies varying from a regular ring to a random 
network. We find that for sufficiently strong coupling strengths the trajectories of this system escape to infinity in 
the regular ring topology, i.e. the amplitude of the oscillations grows explosively without bound. We term this 
phenomenon a blow-up. The surprising central result of this work is the following: the unbounded growth is 
suppressed when some of the regular connections are dynamically randomized. Namely, dynamic random links 
successfully reign in the blow-ups in this complex network, and normal limit cycle behavior ensues for sufficiently 
fast network rewiring. 

This counter-intuitive result is significant in the general context of taming infinities, as infinities arising in 
theoretical constructs of physical phenomena deeply concern scientists across disciplines'". When physical quant- 
ities tend to infinite values in a model, it is often considered an indication of the failure of the theory, or symptoms 
of the limitations of the modeling. Removing infinities, for instance through theoretical tools such as renorma- 
lization, or through improvement of the underlying model of the physical process, has had much impact on 
physics in general. Further from the point of view of potential applications it is indeed crucial that no quantity 
grows without bound. So the drive to find mechanisms that can quell infinities in measurable quantities stems 
from both methodological and pragmatic motivations'". Specifically, our observations in this work may yield 
methods to control and prevent blow-ups in extended oscillator systems which are widely employed in a 
multitude of engineered devices. Lastly, a certain degree of randomness in spatial links may be closer to physical 
reality in many systems of biological, technological and physical significance'. So our finding that random links 
prevent unphysical blow-ups suggests an underlying mechanism by which complex systems can avoid explosive 
growths. 

Results 

We start with a generic coupled system comprised of nonlinear local dynamics and a term modeling the coupling 
interaction. The isolated (uncoupled) dynamics at each node of the lattice/network is given by X = F(X) , where X 
is a m-dimensional vector of dynamical variables and F(X) is the velocity field. The general form of such coupled 
dynamical elements is given by following equation: 



X,- = F(X,)-he£/i,H(X,Xj-; 



i=l,. 



,N, 



(1) 
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Figure 1 | Evolution of x,(t) (i = 1,...100), after a transient time of 15000 (left), and the limit cycles of three representative oscillators stacked one 
above the other (right), for network rewiring period r equal to 0.01 (top row) and 0.1 (bottom row). Here coupling strengthe = 0.6 which is significantly 
greater than critical coupling strength (e^ = O.Sfor these parameter values) and the fraction of random links p is 0.8. Blow-up is prevented in both 
the cases, but faster rewiring leads to a synchronized state (top row) whereas slower rewiring leads to de-synchronization and distortion in the shape of the 
limit cycles (right panel, bottom row). Note that for this value of coupling strength, regular nearest neighbor coupling {p = 0) would lead to a 
blow-up. 



where /,y are the elements of a connectivity matrix, and need not in 
general be symmetric". The coupling strength is given by e and H(X„ 
Xy) is the coupling function determined by the nature of interactions 
between dynamical elements and j. 

At the nodes, we first consider the prototypical Van der Pol oscU- 
lator'''" with nonlinear damping governed by the second-order dif- 
ferential equation: 

x + ^i{x^-l)x + x = 0 (2) 

where x is the dynamical variable and /( is a parameter which deter- 
mines the nature of the dynamics. The relaxation oscillations arising 
in this system, for fi> 0, are characterized by slow asymptotic beha- 
vior and sudden discontinuous jumps. Associating i=y gives: 

k=f{x,y)=y 

y^g{x,y)^-n{x^-'^)y-x 

The Van der Pol equation given above is relevant in modeling phe- 
nomena in physical and biological sciences. Extensions of such oscil- 
lators have been used to model the electrical activity of the heart and 
the action potentials of neurons". The equation has also been utilised 
in seismology to model the two plates in a geological fault". 

Now we consider a ring of N nonlinearly coupled Van-der-Pol 
oscUlators'^"^", namely a specific form of Eqn. 1 with m = 2, X, = 
{Xn yts, F = ^} and H(X,-, X^) = {/(%,, y^ - /{(x„ yd, g(xj, yj) - g{Xi, 
yi)} given as follows: 

Xi =f{xi,yi) + I \f{xi-i ,yi- 1 ) +f{xi+ 1 1 ) - 2f{xi,yi)] 

\ (4) 
h =g{xuyi) + - +^(xi+i,>'i+i) -2^(xi,7i)] 

Here index ; specifies the site/node in the ring, with the nodal on-site 
dynamics being a Van der Pol relaxation oscillator. Starting with this 
regular ring, we construct small world networks^ by rewiring with 



probability p (0 £ p £ 1), the nearest neighbor links of a node to 
some non-local node chosen at random. Additionally the links are 
dynamic, which implies that the specific set of nearest neighbour 
links that get rewired to random nodes varies from time to time. 
So the underlying web of connections changes over time^'"^'', and 
such dynamic rewiring is expected to be widely prevalent in response 
to environmental influences or internal adaptations^^"^'. We imple- 
ment this by rewiring the network at time intervals of r (see meth- 
ods). So a new connectivity matrix having the same proportion of 
random and nearest neighbor connections is formed after every time 
r. So r gives the time-scale of network switches and as r becomes 
larger the network changes become more infrequent, with the limit of 
r ^ 00 corresponding to the standard scenario of quenched links. 

In order to get the dynamics of the complete system, we solved the 
2N coupled non-linear ordinary differential equations (ODE) given 
by Eqns. 3-4 with dynamic network topologies for various coupling 
strengths and system sizes. We studied the shapes and sizes of the 
limit cycles by varying the fraction of random links p, and network 
switching time period r. 

We first describe our principal observations for oscillators coupled 
to two nearest neighbors on a periodic lattice (i.e. p = 0) for various 
coupling strengths. When the coupling strength is between 0 and a 
critical coupling strength e^, the system yields regular bounded beha- 
vior. However when the coupling strength exceeds this value, i.e. e > 
e^, there is a blow-up in the system. Namely, the amplitude of the 
oscillations grows explosively in an unbounded fashion^"'^" For 
instance, in a typical blow-up the amplitude grows from O (1) to 
0(10-*) in time ~ 10"'. 

Now, very different behavior from that described above ensues 
when the nearest neighbor links are dynamically rewired. We show 
the results for representative values of p and r in Figure 1. Notice that 
the coupling strength shown in the figure is greater than the critical 
value Cc, i.e. for this coupling strength the limit cycles grow explo- 
sively for regular nearest neighbor coupling. However, it is clearly 
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evident from the figure that the blow-up has been tamed when some 
fraction of the nearest neighbor Hnks are replaced by random links 
(i.ep > 0). Namely, for sufficiently random networks all oscillators 
remain confined in a bounded region of phase space. 

To quantitatively characterize the bounded-unbounded transi- 
tions, we introduce a "boundedness order parameter" Zho^nd given 
simply as the fraction of initial conditions for which the dynamics 
remains bounded. Thus Zjjoimd gives the probability of successful 
prevention of blow-ups. Fig. 2 displays this quantity over a range 
oip and e. It is evident from Fig. 2 that we always get a bounded state 
if the fraction of random links p is greater than a critical value. 

The fraction of random links necessary to obtain bounded 
dynamics with probability close to 1 depends on the network switch- 
ing time period r and coupling strength e as indicated in Fig. 2. In 
general, for fast switching topologies, namely small r, the critical 
number of random links needed to suppress a blow-up is small. 
However, for very slow network changes one may not be able to quell 
explosive growth even with large fraction of random links. So it is 
clear that rapidly changing links are the key to effectively control the 
dynamics. 

Further, we studied the case of undirected networks (where / in 
Eqn. 1 is a symmetric matrix). Qualitatively similar results are 
obtained though with reduced range for blow-up control (see lower 
panel of Fig. 2). 

To broadly rationalize this behavior we consider the dynamics of 
an individual oscillator in the system as a sum of the nodal dynamics 
and the coupling contributions of the form Ci{t) and C2{t), whose 
strength is given by coupling constant e. 



l + A/l-4^ie2ci(c2-/(Ci) 



x=y + eci[t) 

y= — ^i[x^ — l)y — x-\-ec2{t) 



(5) 



When Ci = 0 and C2 = 0 the equations reduce to a single Van der-Pol 
oscillator, with a stable limit cycle and an unstable spiral repellor at 
the origin (x* = 0, = 0). When Ci 0, Cj ^ 0, the effective 
dynamics of each oscillator in the network (Eqn. 5) can be seen as 
an oscillator with fluctuating parameters Ci and 02^". Note that when 
Ci 7^ 0 one obtains the fixed point solutions: 



ICfiCi 



(6) 



y 



-eci 



This implies that we get two fixed point solutions for each oscillator 
in the coupled system. For very weak coupling (e-^ 0) the additional 
fixed point sits at infinity. As coupling strength increases this fixed 
point moves towards the limit cycle, and at some critical coupling 
strength collides with it, undergoing a saddle node bifurcation. This 
global bifurcation leads to the destruction of the attracting limit cycle 
and consequent blow-up. Now for noisy parameters Ci and C2, such a 
collision leading to annihilation of the limit cycle occurs for suffi- 
ciently large fluctuations in Ci and C2™- Since the magnitude of fluc- 
tuations is governed directly by coupling strength, this implies that 
blow-ups may occur for sufficiently large e. 

In this scenario, the underlying mechanism that controls the 
blow-up is the following: when the links are changed rapidly (i.e. r 
is small) to first approximation the terms Ci{t) and C2(0 in Eqn. 5 
averages to zero, as it is a sum of uncorrected random inputs. 
This effectively reduces Eqn. 5 to that of a single Van der Pol 
oscillator for which there is no blow-up. Alternately, one can 
argue that for rapid random switching of coupling connections 
the deviations from the single oscillator dynamics is essentially a 
sum of uncorrected contributions, which to first approximation 
average out to zero. 

Synchronization of the Bounded State 

Returning to the space-time plots displayed in Fig. 1, we notice an 
additional feature. In the figure, two cases with different network 
rewiring time periods are shown: a fast rewiring case and a slow 
network change scenario. While there is no blow-up in either case, 
the emergent bounded state is different for networks rewiring at 
different rates. When the underlying network changes rapidly a syn- 
chronized state is obtained'^ as evident in the representative 
example in the upper panel of Fig. 1 . However, slow network changes 
yield bounded dynamics that is not synchronized as seen in the 
representative example in the lower panel of Fig. 1. 
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Figure 2 | Transition between bounded and unbounded dynamics in a system of 100 coupled oscillators, as reflected in the variation of the 
boundedness order parameter 2(,„„„^ (see text) with increasing coupling strength e and fraction of random links p, for directed networks with network 
switching time period r equal to (A) 0.01, (B) 0.1 and (C) 1, and for undirected networks with r equal to (D) 0.01, (E) 0.1 and (F) 1. 
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We investigate in detaO the synchronization properties of the 
system. Specifically we attempt to answer the following question: 
when random coupling prevents a blow-up, does it yield a state that 
is synchronized, or not? In order to quantify the degree of synchron- 
ization in the network we compute an average error function as the 
synchronization order parameter, defined as the mean square devi- 
ation of the instantaneous states of the nodes: 

i=l 

where (x{t)) is the space average at time t. This quantity Z, averaged 
over time N and over different initial conditions, is denoted as Zsy„c- 
When 2s^„c = 0 we have complete synchronization. 

Note that the time needed to reach synchronization is not simply 
monotonic with respect to fraction of random links"". Instead we 
observe a region, at intermediate values of p, where the transience 
becomes large. This kind of behavior is reminiscent of the long 
persistence of transience arising for intermediate values of network 
connectivity'^. So in our calculation of the relevant order parameters 
we have left sufficiently long transience ( — 15000). 

It is evident from Fig. 3 that in case of fast changing networks (i.e. r 
is sufficiently small) the bounded state that arises after random rewir- 
ing is also synchronized. This is evident quantitatively through the 
vanishingly small synchronization error (Zsy„c ~ 0) in networks 
where blow-ups have been completely suppressed (i.e. Zhound = !)■ 
On the other hand, for slowly changing networks the synchroniza- 
tion error is significant, even when the dynamics is bounded i.e. Zsy„c 
is not close to zero even when Zi,^^„j = 1, as evident in Fig. 3B. So 
infrequent network changes have led to a state which is bounded but 
not synchronized, i.e. slowly switching random links are effective in 
reigning in the unbounded dynamics, but they do not manage to 
produce synchronization. 

This relaxation oscillator network then yields three kinds of 
dynamical states: (i) bounded synchronized motion; (ii) bounded 
unsynchronized dynamics; and (iii) unbounded dynamics. The nat- 
ure of the dynamics depends on the interplay of the coupling strength 
e, rewiring probability p and network switching time r. 

Probabilistic Switching of Links 

Periodic switching of links occurs in situations where the connec- 
tions are determined by some global external periodic influence. 
However this is not the most realistic scenario, as the interaction 
patterns may not change periodically in time. In such cases we must 



consider a probabilistic model of link switching, such as in^*-'* ''. In 
this subsection, we study such randomly switched networks in order 
to determine if the ability to tame explosive growth is robust to the 
manner in which the links change. Namely, we verify the generality 
of above observations by investigating the dynamics on a network 
whose links switch randomly and asynchronously in time. 

Specifically, a link has a probability of being rewired (the "link 
rewiring probability") in some time interval (see methods). So when 
pr is higher, the links are switched around more often. Note that as 
before, when a node rewires, it links to its nearest neighbors with 
probability 1 ^ p and to some random site with probability p. 

The results are displayed in Figs. 4-6. As in the case of periodically 
changing networks, here too random links manage to control the 
explosive growth in the system for sufficiently large fraction of ran- 
dom links. When the probability of rewiring the links is larger, 
namely when the network changes more often, the critical coupling 
strength at which blow-up occurs is pushed back, yielding a larger 
region of bounded dynamics. 

The interplay of the fraction of random links p and the link rewir- 
ing probability p^ reveals some interesting scaling relations. Notice 
first that the critical coupling strength beyond which blow-ups 
occur, scales with the link rewiring probability p^ and fraction of 
random links p as: 

^c-{pPrY (8) 

Fitting the data using the nonlinear least-squares Marquardt- 
Levenberg fitting algorithm yields /? = 0.1 19 ± 0.001. This is clearly 
evident quantitatively in the data collapse displayed in Fig. 5A. This 
scaling relation implies that as the connections change more fre- 
quently, and there is larger fraction of random links in the connec- 
tivity matrix, the range over which bounded dynamics is obtained 
becomes larger. Also notice that the effective quantity that is relevant 
in the scaling relation is the product ofp andpr> i-e.p andp^ may vary, 
but the emergent phenomena is the same if their product remains the 
same. 

Further calculating the minimum fraction of random links pc 
necessary to prevent blow-ups, for a fixed coupling strength, at dif- 
ferent link rewiring probabilities p^ reveals that 

Pc-(Prf (9) 

The results are displayed in Fig. 5B for two different values of coup- 
ling. Using the nonlinear least-squares Marquardt-Levenberg fitting 
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Figure 3 | Variation of the Boundedness Order Parameter Zh„^„j (solid red curve) and the Synchronization Order Parameter Z^^^ (dotted blue curve), 
averaged over 100 initial conditions for varying fraction of random links p. Here system size JV = 100, coupling strength e = 0.6 and network rewiring 
time period (A) r = 0.01 and (B) r = 0.1. 
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Figure 4 | Transition between bounded and unbounded dynamics in a system of 100 coupled oscillators with probabilistically switched links, as 
reflected in the variation of the boundedness order parameter Zi,„^„j (see text) with increasing couphng strength e and fraction of random links p. The 

link switching probabiUty p, is: (A) 0.1, (B) 0.3 (C) 0.6 and (D) 0.9. 



algorithm one obtains i5 = -1.04 ± 0.008 for e = 0.7 and ^ = -0.98 
± 0.008 for 6 = 0.9. 

So when the system has large spatial randomness infrequent 
changes of links manage to suppress blow-ups, while in a system 
with few spatially random connections the links need to change more 
often in order to contain the blow-ups. This again implies (as did the 
scaling form in Fig. 5 A) that blow-ups can actually be tamed at quite 
small link rewiring probabilities and in networks with small fraction 
of random Unks, underscoring the broader relevance of our observations. 

Lastly, we observe that for networks that change connections fast 
(i.e. when is large) the bounded state that arises after random 



rewiring is also synchronized, as evident quantitatively through the 
vanishingly small synchronization error in these networks (cf 
Fig. 6A). However for slowly changing networks (i.e. when is 
small) the synchronization error is quite significant even when the 
dynamics is bounded (cf Fig. 6B), namely, slowly changing links 
yield a state that is bounded but not synchronized. 

Generality of our central result 

Finally, in order to further gauge the generality of the above findings 
we investigated different coupling forms and local dynamics. For 
instance, we studied some other forms of nonlinear coupling, such as: 
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Figure 5 | (A) Variation of the critical coupling strength at which blow-up occurs, e^, with respect to the product of the link rewiring probability and 
fraction of random links p. The points are data from a wide range of p and p^, and they collapse to the scaling form shown by the line: ec~(p prY with P = 
0.119. (B) Variation of the minimum fraction of random links necessary to prevent blow-ups vs link rewiring probability p^. at coupling strengths e = 0.7 
(red) and 0.9 (blue). The straight lines show the fit, using the nonlinear least-squares Marquardt-Levenberg fitting algorithm, to the scaling form: p^ ' ■ 
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Xi =f (^Xi +^[xi_i+Xi+i-2Xi],yi+^[yi_i+yi+i-2yi]J 



(10) 



and 



^. =/(^My,) + - {x]_ 1 + 1 - 2xf ] ) 

h =g{x,.y:) + \ {y]-i +yhi -27?] ) 



(11) 



Now, simulating all these different coupling forms above (Eqns. 10 
and 11), with some fraction of links randomized, displayed the same 
qualitative phenomena: After a certain critical coupling strength we 
encountered a blow-up in a regular ring of these oscillators, with the 
blow-ups occurring very early in some systems (e.g. e^- < 0. 1 in Eqn. 
11). However as observed in the earlier systems, dynamically ran- 
domizing the links in the network suppressed this explosive growth 
yielding bounded dynamics over a larger range of coupling strengths. 
Note that the range of the bounded region depends on the specifics of 
the coupling form and the local dynamics, as well as the fraction and 
frequency of change of the random links. 



Lastly, we have also investigated the global behavior of the system 
under different nodal dynamics. For instance, we investigated collec- 
tions of nonlinearly coupled Stuart-Landau oscillators, where the 
local dynamics was given by: 



Iz|')z 



(12) 



where co is the frequency and z{t) = x{t) + iy{t). 

We also studied nonlinearly coupled networks of FitzHugh- 
Nagumo oscillators, modeling neuronal populations, described by 
the local dynamics'^-'"': 



x — x~x^/3—y + I 
y = ^{x + a-by) 



(13) 



where x is the membrane potential, y is the recovery variable and / is 
the magnitude of stimulus current. 

Further, we have explored the heterogeneous case where the non- 
linearity parameter ji is distributed over a range of positive values for 
the local Van der Pol oscillators. 

Our extensive numerical simulations showed qualitatively similar 
behavior for all of the above systems (see Fig. 7 for representative 
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Figure 7 | (A) Phase portraits and (B) time evolution of the membrane potential, for two representative FitzHugh-Nagumo model neurons in the 
network given by Eqn. 13, with fraction of random links p = 0.5, time period of changing the links r = 0.01, and coupling strength e = 0.54 (magenta and 
blue curves). The case of a static regular network, exhibiting blow-up, is shown alongside as well (red curve). Note that the bounded orbits obtained under 
dynamic random links are not synchronized (magenta and blue curves). 
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examples). This strongly suggests that the phenomenon of preven- 
tion of blow-ups through dynamic random links is quite general. 

Conclusions 

In summary, we have studied the dynamics of a collection of relaxa- 
tion oscillators under varying coupling topologies, ranging from a 
regular ring to a random network. We find that the system blows up 
under regular ring topology, for sufficiently strong coupling strength. 
However when some fraction of the links are randomized dynam- 
ically, this blow-up is suppressed and the system remains bounded. 
So our results suggest an underlying mechanism by which complex 
systems can avoid a catastrophic blow-up. Further from the stand 
point of potential applications, our observations indicate a method to 
control and prevent blow-ups in coupled oscillators that are com- 
monplace in a variety of engineered systems. 

Methods 

We present the method for changing the connectivity of the network. 
Algorithm for periodically switching networks 

1. We start with a 1 dimensional ring with iV nodes. Each node i (i — 1, ... JV) has 2 
Unks corresponding to its two nearest neighbors. These regular links of the node 
are cut with probability p and connected to some randomly selected distant 
node. In the network, on average, a fraction p of the links are to random nodes 
and fraction (1 — p) are to nearest neighbors. So this network corresponds to a 
standard small world network^ with degree 2 and fraction of random links given 
by p. 

2. The above process of generating a small world network of degree 2 and fractionp 
of random links is repeated after time interval r, i.e. the network changes to this 
new configuration after time r. When the connections switch, some nearest 
neighbor links of the previous network configuration may change to random 
links and some random links may switch back to regular ones. 

So the network switches from one configuration to another periodically, keeping the 
qualitative nature of the network invariant. Namely, at any instant of time the con- 
nectivity matrix has a fraction (1 — p) of regular links and fraction p of random links. 

Algorithm for stochastically switching networks 

1 . We again start with a 1 dimensional ring with N nodes, where each node f { f — 1 , 
... N) has 2 links. Again fractionp of the links are to random nodes and fraction 
(1 — p) are to nearest neighbors. 

2. Now each individual node i has probability p^ of changing its links in a specified 
time interval t. 

3. If the links of a certain node are selected for rewiring, they rewire such that with 
probability p they connect to random neighbors and with probability ( 1 — p) 
they connect to nearest neighbors. 

In this algorithm, the nodes change their links independently and stochastically. The 
probability of the new links being random or regular is determined by p as in the first 
algorithm (and as in the standard small-world scenario). Note that in the first algo- 
rithm r is the time period of link changes which occur globally throughout the 
network, while in the second algorithm uncorrelated changes occur at the local level. 

Our simulations show that these two very different ways of varying the network 
both yield the same qualitative phenomena: namely changing links helps to suppress 
blow-ups in the system. 

Method for simulating the dynamics of the system. We solved the 2N coupled non- 
linear ordinary differential equations (ODE), given by Eqns. 4, 10, 1 1 by integrating 
the ODEs with the fourth order Runge-Kutta algorithm. We have checked that the 
results are qualitatively same over a large range of step sizes (—10"' to 10"^), in order 
to rule out artefacts arising from the stiffness problem in such differential equations. 

We obtained the shapes and sizes of the limit cycles arising in this network, for 
system sizes ranging from N — 10 to N — 10^, under varying fraction of random links 
p (0 <p < 1). We investigated a large range of network switching time periods r (0.01 
< r < 1) for the case of periodically switched networks, and 0 < p^ 1 for the case of 
stochastically switched networks (with t — 0.001). 

System parameters. In the system given by Eqns. 3 and 4 we have investigated a range 
of nonlinearity parameters f.i and we present the results for the representative case of fi 
= 1.5 in the examples displayed in this paper. In the system given by Eqn. 13, 
parameters 0, a, b are taken to be 0.08, 0.7, 0.8 respectively. 
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